Comprehensive bioinformatic analysis reveals a cancer-associated fibroblast gene signature as a poor prognostic factor and potential therapeutic target in gastric cancer

Background Gastric cancer is one of the deadliest cancers, currently available therapies have limited success. Cancer-associated fibroblasts (CAFs) are pivotal cells in the stroma of gastric tumors posing a great risk for progression and chemoresistance. The poor prognostic signature for CAFs is not clear in gastric cancer, and drugs that target CAFs are lacking in the clinic. In this study, we aim to identify a poor prognostic gene signature for CAFs, targeting which may increase the therapeutic success in gastric cancer. Methods We analyzed four GEO datasets with a network-based approach and validated key CAF markers in The Cancer Genome Atlas (TCGA) and The Asian Cancer Research Group (ACRG) cohorts. We implemented stepwise multivariate Cox regression guided by a pan-cancer analysis in TCGA to identify a poor prognostic gene signature for CAF infiltration in gastric cancer. Lastly, we conducted a database search for drugs targeting the signature genes. Results Our study revealed the COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC as the key CAF markers in gastric cancer. Analysis of the TCGA and ACRG cohorts validated their upregulation and poor prognostic significance. The stepwise multivariate Cox regression elucidated COL1A1 and COL5A1, together with ITGA4, Emilin1, and TSPAN9 as poor prognostic signature genes for CAF infiltration. The search on drug databases revealed collagenase clostridium histolyticum, ocriplasmin, halofuginone, natalizumab, firategrast, and BIO-1211 as the potential drugs for further investigation. Conclusions Our study demonstrated the central role of extracellular matrix components secreted and remodeled by CAFs in gastric cancer. The gene signature we identified in this study carries high potential as a predictive tool for poor prognosis in gastric cancer patients. Elucidating the mechanisms by which the signature genes contribute to poor patient outcomes can lead to the discovery of more potent molecular-targeted agents and increase the therapeutic success in gastric cancer. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-022-09736-5.


Background
Gastric cancer is the fifth most common cancer worldwide and the fourth leading cause of cancerrelated deaths, the GLOBOCAN 2020 statistics report. More than one million people were diagnosed with gastric cancer, and more than 750,000 deaths occurred due to gastric cancer in 2020 [1]. Stomach adenocarcinomas (STAD) constitute almost 95% of all gastric cancer cases. The mainstay of treatment in localized stomach adenocarcinoma is gastrectomy with total lymphadenectomy and chemotherapy [2]. However, the tumors are commonly metastatic at the time of diagnosis. At this stage, complete resection is impossible, and currently available chemotherapeutics fail due to chemoresistance [2,3].
Molecular-targeted agents against tumor-specific biomarkers and immunotherapy increased the treatment efficacy in certain cancers such as breast cancer, lung cancer, and melanoma [4]. However, similar success is not achieved in gastric cancer yet. Currently, molecular-targeted agents, nivolumab (anti-PD-1), pembrolizumab (anti-PD-1), ramucirumab (anti-VEGFR2), and trastuzumab (anti-HER2), are approved in gastric cancer treatment. Unfortunately, they have limited efficacy on the overall survival of a limited group of advanced-stage patients with target positivity [5].
The tumor microenvironment is a great challenge for the treatment of cancer. Dynamic interactions with the extracellular matrix (ECM) and cellular components in the tumor microenvironment potentiate the aggressiveness of cancer cells and limit their response to anti-cancer agents [6]. Cancer-associated fibroblasts (CAFs) are critical components in the cellular compartment of the tumor microenvironment that assemble and remodel the ECM. They originate from activated fibroblasts and the endothelial or epithelial cells undergoing epithelial-mesenchymal transition (EMT). They secrete various ECM proteins and soluble mediators that potentiate pro-tumorigenic signaling pathways via activation of transmembrane receptors -mainly integrins. Moreover, CAFs remodel the ECM to form a protective barrier against immune surveillance and the diffusion of anti-cancer agents [7]. CAF infiltration is associated with a dismal prognosis in gastric cancer [8,9]. CAFs were identified as the greatest risk factor in tumor microenvironment phenotype with the poorest overall survival in gastric cancer patients [10]. Therefore, addressing CAFs is essential for the treatment of gastric cancer. However, the poor prognostic markers for CAF infiltration in gastric cancer are not clear and drugs that target CAFmediated processes are lacking in the clinic.
In this study, we aim to identify the poor prognostic gene signature for CAF infiltration targeting which may increase the therapeutic efficacy in a large group of gastric cancer patients. We comprehensively analyzed four gastric cancer GEO gene expression datasets using a network-based approach and identified key markers for CAF infiltration. After validation of the markers in The Cancer Genome Atlas (TCGA) and Asian Cancer Research Group (ACRG) cohorts, we elucidated a poor prognostic gene signature for CAF infiltration in gastric cancer using stepwise multivariate Cox regression. Lastly, we compiled the list of currently available drugs that may have a therapeutic potential in gastric cancer by targeting the signature genes we identified ( Fig. 1 summarizes the steps followed in the study).

Data collection and identification of differentially expressed genes in gastric cancer
We analyzed four GEO expression profiling datasets (GSE13911, GSE29272, GSE79973, GSE118916) [11][12][13][14], which used Affymetrix Human Genome or Gene Expression Arrays for profiling (https:// www. ncbi. nlm. nih. gov/ geo/). All four datasets bear a comparable number of samples from gastric cancer tissues and noncancerous gastric tissues. Additional file 1: Table S1 lists the number of tissue samples and the profiling platforms in each dataset. In total, we analyzed 200 non-tumor gastric and 207 gastric tumor samples. To identify differentially expressed genes (DEGs) in gastric tumors compared to non-cancerous stomach samples, we used the GEO2R web tool (https:// www. ncbi. nlm. nih. gov/ geo/ geo2r/). We applied log transformation and Benjamini & Hochberg's (False discovery rate) method to adjust the p-values (p-value significance cut-off = 0.01). The genes were filtered based on their log2-fold change (logFC) values. We accepted the genes with the log FC value >1 as the upregulated genes and with the log FC value < −1 as the downregulated genes. After identifying the DEGs in each dataset, we performed Venn Analysis to find DEGs common to all four datasets using the jvenn (an Ucaryilmaz Metin and Ozcan BMC Cancer (2022) 22:692 interactive Venn diagram viewer) (http:// jvenn. toulo use. inra. fr/ app/ index. html).

Functional annotation and enrichment analysis
To perform functional enrichment and annotation clustering analysis of the DEGs, we used The Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Version 6.8) [15] (https:// david. ncifc rf. gov/). To understand the cellular compartments (GO-CC), molecular functions (GO-MF), and biological processes (GO-BP) at which the DEGs enriched, we performed gene ontology (GO) analysis. To understand the pathways at which the DEGs operate, we performed the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis [16]. The cut-off value for significance was chosen as p < 0.05. We analyzed gene lists for upregulated genes, downregulated genes, and all DEGs separately.
To compare the gene identities and gene ontologies enriched in different lists, we used Metascape [17] (https:// metas cape. org). To dissect the similarities and dissimilarities between gene lists, we analyzed the Circos plots, clustering dendrograms, network layouts for enriched gene ontologies, the proportion of the genes from different lists that fall into the same gene ontologies, and enrichment p-values.

Protein-protein interaction network analysis
To identify the central CAF markers and assess their potential connections and interactions with the protein products of other DEGs in gastric cancer, we constructed the protein-protein interaction (PPI) network of the DEGs using The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING Version 11.0) [18] (https:// string-db. org/). We utilized the minimum required interaction score as high confidence (0.7). We analyzed the resulting PPI network in Cytoscape (Version 3.8.2) to infer the topological parameters of each node [19] (https:// cytos cape. org/). To investigate the network modules, we used Molecular Complex Detection (MCODE) plugin at Cytoscape. To construct a local network for key CAF markers and their first neighbors, we used the Cytohubba plugin at Cytoscape [20,21]. To investigate the interactors of ITGA4 we searched inBio Discover ™ by Intomics A/S (https:// inbio-disco ver. com/) (Intomics A/S has not endorsed the results of the published article) [22].

Gene expression profiling
To confirm the differential expression of the CAF markers in gastric cancer, we comparatively analyzed the gene expression profiles of 34 non-cancerous gastric tissues and 415 stomach adenocarcinoma samples in the TCGA dataset using the UALCAN (University of Alabama Cancer Database) (http:// ualcan. path. uab. edu/) [23]. We also investigated the differential expression of the CAF markers by tumor grade and stage (the unpaired t-test was used for statistical analysis). To validate the results from the UALCAN, we analyzed The Genotype-Tissue Expression Project (GTEx) data on GEPIA2 [24] (http:// gepia2. cancer-pku. cn).
To investigate the differential expression of CAF markers in diffuse vs. intestinal subtypes and mesenchymal vs. epithelial phenotypes of gastric adenocarcinoma, we analyzed the ACRG cohort on GEO2R (GSE66229) [25]. We used the same set of parameters to identify DEGs in GSE66229 as for the datasets GSE13911, GSE29272, GSE79973, and GSE118916. Then we analyzed the expression profile graphs of all patients for the six CAF markers. To investigate the differential expression of CAF markers in other cancers we analyzed the TCGA data on Tumor Immune Estimation Resource 2.0 (TIMER2.0) [26] (https:// timer. cistr ome. org). Then, we extracted the pancancer expression profile graphs for CAF markers.

Survival analysis
To understand the impact of the CAF markers on the survival of gastric cancer patients, we performed the Kaplan-Meier (KM) survival analysis of TCGA stomach adenocarcinoma samples on TIMER 2.0. The stomach adenocarcinoma samples split into high or low expression groups based on the median expression level for each gene. Additionally, KM-Survival Curve for COL1A2 was extracted from the UALCAN to assess its prognostic role in gastric cancer. We generated the heatmaps that show the z-scores for each gene in distinct cancers using the "gene outcome" module in TIMER 2.0. We constituted the KM-survival curves for CAF infiltration that integrate gene expression data with the "immune association" tool in TIMER 2.0, which utilizes the log-rank test. Then we extracted the hazard ratios (HR), z-scores, and p-values from the multivariate Cox proportional hazard regression models built on TIMER 2.0.

Gene correlation analysis
We investigated the correlation between the individual gene expression and CAF infiltration in different cancers by the "immune association" tool in TIMER 2.0. To examine the correlation between distinct genes, we used the "Gene Correlation" module in TIMER 2.0, which gives purity adjusted correlation coefficients calculated by partial spearman rank correlation. We extracted the correlation coefficients and the p-values from TIMER2.0.

Data visualization
To draw the bubble plots for enriched ontology terms, we used the BIC imageGP tool. To plot the gene expression profiles and hazards ratios for CAF infiltration and gene expression in different cancers, we used GraphPad Prism9.

Identification of cancer-associated fibroblast markers in gastric cancer
To identify the CAF markers in gastric cancer we first investigated the DEGs in gastric cancer vs. noncancerous gastric tissues. We analyzed GSE13911, GSE29272, GSE79973, and GSE118916 datasets which included microarray data for gastric tumors of various subtypes. To determine the most relevant markers, we applied strict criteria for the identification of the DEGs. We accepted genes with a log FC value of >1 or < −1 and an adjusted p-value <0.01, rather than accepting all genes with a log FC different than zero and p-value <0.05 as DEGs. Figure Table S2).
To illuminate the biological functions and the pathways the DEGs enrich, we performed the functional annotation and enrichment analysis of 83 DEGs. The functional annotation clustering with the highest classification stringency in DAVID revealed 3 clusters ( Table 1). The cluster with the highest enrichment score included the ECM-receptor interaction, focal adhesion, and PI3K-Akt signaling pathway.
Functional enrichment analysis of the upregulated genes indicated a key role in ECM organization, ECM remodeling, ECM-receptor interaction, and activation of pro-tumorigenic signaling pathways (Fig. 2e, Additional file 1: Table S3). The most enriched molecular functions were ECM structural constituent, platelet-derived growth factor binding, and integrin-binding. The top KEGG pathways related to the upregulated genes were ECM-receptor interaction, focal adhesion, and PI3K-Akt signaling (Fig. 2g, Additional file 1: Table S3). The downregulated genes enriched in metabolic processes and ion homeostasis ( Fig. 2f-g, Additional file 1: Table S4). These findings pointed out the upregulation of the ECM organization and remodeling in gastric cancer and suggested the involvement of CAFs.
To determine the CAF markers upregulated in gastric cancer we comparatively analyzed the list of upregulated DEGs in gastric cancer with an extended list of CAF markers, in terms of identity and gene ontology in

Protein-protein interaction network analysis and identification of the key CAF markers
To identify the key CAF markers and investigate their interactions with the protein products of DEGs in gastric cancer, we constructed a PPI network in STRING (Fig. 3a). The analysis of this network on Cytoscape 3.8.2. revealed a prominent hub composed almost totally of ECM components and CAF markers. The top six proteins with the highest degree in this hub were all CAF markers: COL3A1 (collagen type III alpha 1 chain), FN1 (fibronectin 1), COL1A2 (collagen type I alpha 2 chain), COL1A1 (collagen type I alpha 1 chain), COL5A1 (collagen type V alpha 1 chain), and SPARC (cysteinerich acidic matrix-associated protein) respectively (Fig. 3b). Table 2 shows the topological parameters for these markers. The topological parameters for the whole PPI network are listed in Additional file 1: Table S5. The components of the hub were all upregulated genes in gastric tumors, except for COL2A1 (collagen type II alpha 1 chain) (Additional file 2: Fig. S1). These findings strengthened the central role of CAF markers in gastric cancer.

Analysis of the network modules CAF markers involved
To identify the modules that the six key CAF markers function, we used the MCODE tool in Cytoscape. MCODE identified three modules. The upregulated ECM protein hub was represented by modules 1 and 2 in Fig. 3c. Five out of 6 CAF markers, COL1A1, COL1A2, COL3A1, COL5A1, and SPARC, were components of module 1, while the FN1 was in module 2.
Then we performed functional enrichment analysis to identify enriched KEGG pathways at each module. "ECM-receptor interaction" and "focal adhesion" were the common enriched KEGG pathways in modules 1 and 2 (Additional file 1: Table S6). "PI3K-Akt signaling pathway" was the third enriched pathway in module 1. Hence the module analysis strengthened the connection of the six CAF markers: COL1A1, COL1A2, COL3A1, COL5A1, FN1 and, SPARC , with the 3 KEGG pathways: ECM-receptor interaction, focal adhesion and, PI3K-Akt signaling pathway in gastric cancer. COL1A1, COL1A2, COL3A1, COL5A1, and FN1 are abundant structural proteins at the ECM. COL1A1 and COL1A2 are produced mainly by fibroblasts and together constitute the type I collagen in the connective tissue. COL3A1 and COL5A1 are the alpha-1 chains of type III and V collagen, which are found in connective tissue together with type I collagen [28,29]. FN1 is a glycoprotein involved in cell adhesion, wound healing, and metastasis. Besides their structural role, these proteins bind to the integrins on the cell membrane, and through focal adhesion kinases, they activate intracellular signaling pathways such as PI3K-Akt and MAPK pathways [30]. SPARC encodes the cysteine-rich acidic matrix-associated protein that is an essential protein for ECM remodeling. It binds to collagens and fibronectin; and regulates the interactions of cells with the ECM [31].
We analyzed the first neighbors of these CAF markers in our PPI network using the Cytohubba tool in Cytoscape. All these CAF markers highly interacted with other structural ECM components: BGN (biglycan), THBS1/2 (thrombospondin 1/2), and VCAN (versican), or ECM remodeling enzymes like SERPINH1 (serpin Table 1 Functional annotation clustering of 83 differentially expressed genes in gastric cancer. Analysis was performed with the highest classification stringency in DAVID (ECM: Extracellular Matrix, GOTERM: Gene Ontology Term, GOTERM BP: GO-biological process, GOTERM MF: GO-molecular function, and GOTERM CC: GO-cellular compartment, KEGG Pathway: Pathways listed in Kyoto Encyclopedia of Genes and Genomes)  (Line colors indicate the type of interaction evidence; edges indicate functional and physical protein associations. Minimum required interaction score: high confidence (0.7). Disconnected nodes are hidden in the network). B Degree sorted circular network graph constructed in Cytoscape. Black arrows show the two genes with the highest degree in the network. The degrees of the nodes decrease counterclockwise. C Modules in the protein-protein interaction network of 83 differentially expressed genes in gastric cancer. Module analysis was performed using the MCODE tool in Cytoscape. The scores, number of nodes, and edges for each module are listed in the D Neighbors of the six CAF markers COL1A1 (orange), COL1A2 (orange), COL3A1 (red), COL5A1 (yellow), FN1(red), and SPARC (yellow). The network was constructed using the Cytohubba tool in Cytoscape (Color grade from red to yellow indicates the descending order of node degrees in the network) family H member 1), SULF1 (sulfatase 1), and TIMP1 (tissue inhibitor matrix metalloproteinase 1) (Fig. 3d).

The correlation of the six CAF markers with CAF infiltration and other CAF markers in gastric cancer
To validate the six key CAF markers we detected with a network-based analysis, we investigated their correlation with CAF infiltration and other CAF markers in gastric cancer. Correlation analysis in TIMER2.0 revealed a statistically significant correlation for all six CAF markers with CAF infiltration in the TCGA stomach adenocarcinoma dataset ( Fig. 4a-f ). The correlation of these markers with the CAF infiltration was even higher than that of poor prognostic CAF signature genes recently identified in gastric cancer: THBS1, THBS2, INHBA (inhibin A), CXCL12 (C-X-C motif chemokine ligand 12), TGFB2 (transforming growth factor-beta 2), VEGFB (vascular endothelial growth factor B), COL10A1 (collagen type X alpha 1 chain), AREG (amphiregulin) and EFNA5 (ephrin A5) (Additional file 2: Fig. S2) [8,32,33]. These findings validated the central role of COL1A1, COL1A2, COL3A1, COL5A1, and FN1 as CAF markers in gastric cancer.
Then we investigated the correlation of the six key CAF markers with the total list of CAF markers at the gene expression level in TCGA stomach adenocarcinoma samples (Fig. 4g). We also added CXCL12, INHBA, THBS1, and THBS2 to the correlation analysis since they are suggested as CAF markers associated with an aggressive phenotype in gastric cancer [32][33][34]. In hierarchical clustering analysis of the correlation matrix, the six CAF markers were highly correlated and clustered with COL11A1 (collagen type XI alpha1 chain), FAP (fibroblast activation protein), INHBA, MMP11 (matrix metalloproteinase 11), S100A4 (S100 calcium-binding protein A4) and THBS2 (Fig. 4g).

Verifying the differential expression and prognostic impact of the six CAF markers in gastric cancer
To validate the differential expression of the six CAF markers in gastric cancer, we analyzed the TCGA data in UALCAN. The expression of all six markers was significantly higher in stomach adenocarcinoma samples, compatible with our results (Fig. 5a). The significant upregulation of these markers in gastric cancer was also verified on GEPIA2 using gastric cancer data from the GTEx dataset (data not shown).
To investigate whether the six CAF markers are involved in gastric cancer progression, we analyzed their differential expression by tumor stage and grade on UALCAN using the TCGA data. For COL1A2, COL3A1, COL5A1, FN1, and SPARC , there was no significant upregulation in stage 1 patients compared to healthy controls (Fig. 5b). However, their expression was significantly higher in stage 2, 3, and 4 samples than in the stage 1 samples. Only for COL1A1, the expression was high starting from stage 1. After stage 2, the increase in COL1A1 expression became much more prominent. These findings suggest that COL1A1 may be involved in both tumorigenesis and tumor progression. On the other hand, COL1A2, COL3A1, COL5A1, FN1, and SPARC may be more involved in tumor progression from stage 1 to stage 2, at which the tumor cells gain the ability to invade surrounding tissues.
The expression of COL1A1, COL1A2, COL5A1, and SPARC was significantly high starting from grade 1  (Fig. 5c). The medians of expression for these markers gradually increased in grade 2 and 3 samples. For COL3A1 and FN1, the upregulation in grade 1 disease compared to healthy gastric tissue was not statistically significant due to the high variance in grade 1. However, their expression was significantly higher in grade 2 and grade 3 compared to grade 1 disease and normal tissue. These findings suggest that all six CAF markers may be associated with a poorly differentiated phenotype in gastric cancer. Then we investigated the impact of upregulated COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC on patient survival with KM-survival analysis (Fig. 5d). High expression of COL1A1, COL3A1, COL5A1, FN1, and SPARC was significantly associated with poor survival in gastric cancer patients (p < 0.05). Although the survival curves for COL1A2 high vs. low expression samples were different, the difference was not significant enough to suggest that the high expression of COL1A2 is associated with poor survival (p = 0.0618). Despite that, the KM-survival curves for COL1A2 high vs. low/ medium expression samples from the TCGA database were statistically different (p = 0.029) (Additional file 2: Fig. S3). These findings support that COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC are associated with poor prognosis in gastric cancer.

The differential expression of the six CAF markers in gastric cancer subtypes
Gastric adenocarcinoma is a heterogeneous disease with two main histopathological subtypes: the intestinal-type and diffuse-type gastric adenocarcinoma. The intestinal type is characterized by organized glandular structures and responds better to chemotherapy. On the other hand, diffuse type is characterized by an undifferentiated phenotype and its response to chemotherapy is dismal [35]. Recent studies for the molecular characterization of gastric tumors also demonstrated that gastric tumors with a mesenchymal phenotype have a worse prognosis and poor response to chemotherapy compared to that with an epithelial phenotype [25]. Since CAF infiltration is associated with EMT in several cancers and CAFs can develop from mesenchymal cells [36][37][38], we asked whether the six key CAF markers are more dominant in gastric tumors with a mesenchymal phenotype. Our analysis in the ACGR cohort revealed that expression of all the six CAF markers was higher in gastric tumors with a mesenchymal phenotype compared with that of epithelial phenotype (Fig. 6). However, the expression profile of these markers did not significantly differ between the diffuse vs. intestinal subtypes (Additional file 2: Fig. S4). Although diffuse-type gastric adenocarcinoma is more associated with a mesenchymal phenotype compared with the intestinal-type [35], it exhibits inter-patient variability in terms of the mesenchymal markers [25]. Therefore, mesenchymal phenotype seems to be a better indicator for the expression of six CAF markers, hence the CAF infiltration.

Pan-cancer analysis of the six CAF markers at the TCGA database
CAFs are the predominant cellular components in the microenvironment of various tumors, especially with a high stroma [39]. To understand whether the six CAF markers we identified in gastric cancer are also upregulated in other cancers we performed a pancancer analysis of these markers in the TCGA dataset. Although each of the markers was upregulated in more than half of the cancers, there were five cancer types besides stomach adenocarcinoma in which all the six markers were upregulated: colon adenocarcinoma (COAD), glioblastoma (GBM), head and neck squamous carcinoma (HNSC), kidney renal cell carcinoma (KIRC) and thyroid carcinoma (THCA) (Fig. 7). Despite that, high expression of all the six markers was not associated with poor prognosis in these cancers (Fig. 8a). CAF infiltration was not even associated with a poor prognosis in colon adenocarcinoma and head and neck squamous carcinoma among these cancers in the TCGA dataset (Fig. 8b).
CAFs display a heterogenous gene expression profile in the stroma of distinct cancers. Hence, they play antitumor or tumor-promoting roles depending on the tumor type [7]. Although the six CAF markers were upregulated in colon adenocarcinoma, glioblastoma, head and neck squamous carcinoma, kidney renal cell carcinoma, and thyroid carcinoma, the overall CAF secretome or inner cellular machinery that responds to the six CAF markers may not end up with a poor prognostic response in these tumors. Therefore, the six CAF markers we identified in gastric cancer may not necessarily indicate a poor prognostic impact on CAF infiltration in these tumors. With further analysis, we identified four other cancers in which all six markers and CAF infiltration were associated with poor prognosis: adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), kidney renal papillary cell carcinoma (KIRP), and mesothelioma (MESO) (Fig. 8a). We comparatively analyzed these cancers with gastric cancer to establish a poor prognostic CAF gene signature with high specificity to gastric cancer.

Prognostic impact of the six CAF markers on gastric tumors with CAF infiltration
CAF infiltration by itself is associated with a hazard ratio of 5.24 in stomach adenocarcinoma in the KM-survival analysis of TCGA data (Table 3). Moreover, the outcome of CAF infiltration worsens with the increasing tumor stage in gastric cancer (Additional file 2: Fig. S5). We investigated whether high expression of the six CAF markers further worsens the prognosis in gastric tumors with CAF infiltration. We performed Cox proportional hazard regression analysis that considers both gene expression profiles and CAF infiltration in the TCGA stomach adenocarcinoma samples. Out of the six CAF markers, the expression of COL1A1, COL1A2, COL3A1, or COL5A1 increased the z-score and hazards ratio for CAF infiltration ( Table 3). The high COL5A1 expression led to the highest increase in the risk for poor survival (z = 2.666, HR = 8.584). The high FN1 or SPARC expression did not increase the z-score and hazards ratio for CAF infiltration.
After that, we investigated whether the high expression of dual combinations of COL1A1, COL1A2, COL3A1, or COL5A1 exacerbates the outcome of CAF infiltration in stomach adenocarcinoma (Table 3). Concomitant high expression of COL1A1 and COL5A1 increased the hazard ratio most for CAF infiltration in TCGA samples (z = 2.924, HR = 11.654). The hazard ratio of CAF infiltration with this dual gene combination was even higher than that with the quadruple combination of collagen subunits. We also investigated the impact of CAF markers highly clustered with COL1A1 and COL5A1 in correlation analysis, namely THBS2, FAP, INHBA, S100A4, COL11A1, or MMP11, on the outcome of CAF infiltration in stomach adenocarcinoma. Except for MMP11, all slightly increased the hazard ratio for CAF infiltration (Additional file 1:  [8,32]. To compare the prognostic significance of these two gene signatures with that of COL1A1 and COL5A1, we investigated the Cox regression models for these two signatures. The z-scores and hazard ratios for both signatures were lower than those for COL1A1 and COL5A1 (Additional file 1: Table S8). These findings suggested a high potential for COL1A1 and COL5A1 as a poor prognostic signature in CAF infiltrated gastric tumors.

Prognostic significance of COL1A1 and COL5A1 for CAF infiltration in other cancers
At the next step, we asked whether COL1A1 and COL5A1 increase the poor prognostic impact of CAF infiltration in four cancers (adrenocortical carcinoma, bladder urothelial carcinoma, kidney renal papillary cell carcinoma, and mesothelioma) with a poor outcome profile for the six CAF markers and CAF infiltration like gastric cancer. Interestingly, the addition of COL1A1 and COL5A1 to the CAF Cox model led to a decreased risk score for CAF infiltration in adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma (Fig. 8c).
Identifying the players for the opposing roles of COL1A1 and COL5A1 in different cancers could reveal new insights into the field. To predict possible players, we extracted the list of genes that correlate with the expression of COL1A1 and COL5A1 in adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma. We identified 25 genes that are highly correlated (rho ≥ 0.5) with both COL1A1 and COL5A1 in all three cancers (Fig. 8d). Then we comparatively analyzed this list with the list of genes upregulated in gastric cancer. The two lists shared seven genes (Fig. 8e, upper circos  plot). The 18 genes common to the three cancers fell into the same gene ontology as the 27 genes upregulated in gastric cancer (Fig. 8e, lower circos plot). Despite these overlaps, more than ten ontologies are differentially enriched in the list of upregulated genes in gastric cancer (Fig. 8f). Among these, "integrin α 4 β 1 pathway" and "peptide crosslinking" were striking, since they were the two ontologies that are also enriched at the upregulated gene list in gastric cancer (See figure on next page.) Fig. 6 The differential expression of six CAF markers in gastric cancers with mesenchymal vs. epithelial phenotype. The differential expression of COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC in gastric cancers with mesenchymal vs. epithelial phenotype and normal gastric tissues from corresponding patients (abbreviated as N-MP: "Normal tissue-MP" for patients with mesenchymal phenotype gastric cancer and N-EP: "Normal tissue-EP" for patients with epithelial phenotype gastric cancer) in the Asian Cancer Research Group gastric cancer dataset (GSE66229). Analysis was performed on GEO2R Ucaryilmaz Metin and Ozcan BMC Cancer compared to the CAF markers list (Fig. 9a). Network layouts for enriched ontology clusters given in Fig. 9b-d better demonstrate that the most striking difference for upregulated genes in gastric cancer was the enrichment of the "integrin α 4 β 1 pathway" compared with the CAF markers list. These observations emphasized the role of the "integrin α 4 β 1 pathway" together with CAF infiltration in gastric cancer and proposed the integrin α 4 β 1 signaling as a poor prognostic factor for CAF infiltration specifically in gastric cancer.

Contribution of Integrin α4β1 pathway to the poor prognostic impact of COL1A1, COL5A1, and CAF infiltration in stomach adenocarcinoma
Integrins are heterodimeric transmembrane proteins that are involved in cell-cell or cell-ECM adhesions. They bind ECM components, mainly collagen, and fibronectin, activate intracellular signaling pathways, and regulate cell survival, proliferation, migration, and differentiation. Integrin α 4 β 1 is a heterodimer of integrin α 4 (ITGA4) and integrin β 1 (ITGB1). ITGB1 couples with a large variety of α integrin subunits [40]. However, ITGA4 couples with integrin β 1 or β 7 subunits [41,42]. The integrin α 4 β 1 , also known as very late antigen-4 (VLA-4), is expressed on various immune cells, mediating the migration of leukocytes to the inflammatory sites via interaction with VCAM-1 (vascular cell adhesion protein 1) [41]. Additionally, it binds to ECM components and takes part in fibronectin assembly [43]. Increased expression of α 4 β 1 integrin is associated with tumor progression and chemoresistance in cancer. The interaction of integrin α 4 β 1 on the tumor cell membrane with the VCAM-1 on vascular endothelial cells is involved in metastasis [41]. The interaction of α 4 β 1 integrin with fibronectin suppressed apoptosis via FAK-mediated suppression of p53, and PI3K/Akt mediated upregulation of Bcl-2 in myeloma cells. Increased integrin α 4 β 1 expression was associated with increased binding of melanoma cells to collagen I and collagen IV, and invasion through fibronectin [44]. Moreover, α 4 integrin was suggested to affect a drug efflux mechanism independent of its coupling with β 1 integrin [45]. However, the mechanisms by which integrin α 4 β 1 heterodimer or α 4 integrin monomer contribute to invasion, metastasis, and chemoresistance in cancer are not exactly known.
To understand whether the integrin α 4 β 1 potentiates the poor prognostic impact of CAFs, we analyzed the Cox regression model for CAF infiltration that considers the expression of ITGA4, ITGB1, or both in addition to COL1A1 and COL5A1. The addition of ITGA4 to the model increased the hazard ratio and z-score in stomach adenocarcinoma (z = 2.963, HR = 12.247) ( Table 4). However, only ITGB1 or ITGA4 and ITGB1 slightly decreased the hazard ratio and z-score, which may be due to a less selective coupling of ITGB1 with several integrin α subtypes (Additional file 1: Table S9). These findings supported that ITGA4 may worsen the poor prognostic impact of CAFs in gastric cancer.
The EMILIN-1 is a member of the elastin microfibrillar interface proteins (EMILINs) family, expressed as a homotrimer at the ECM [50]. Since fibroblasts are the major sources of EMILIN-1 at the ECM, it is accepted as a fibroblast marker [51]. The interaction of EMILIN-1 with integrin α 4 β 1 is involved in cell adhesion and migration [52]. The increase in the poor prognostic impact of CAFs in stomach adenocarcinoma with the addition of EMILIN1 as a covariate to the Cox model was surprising since EMILIN-1 is known as a tumor suppressor that exerts anti-proliferative action via integrin α 4 β 1 in cell and in vivo models [50,53]. Despite that, some reports suggest a pro-tumorigenic role for EMILIN-1 in ovarian serous tumors and osteosarcoma [54,55]. A recent study reported that the action of EMILIN-1 to inhibit the MAPK pathway and suppress proliferation in gastric cancer cells might depend on Tetraspanin9 (TSPAN9) [56], which is a member of the tetraspanin family membrane receptors with four transmembrane domains. These receptors are involved in signal transduction, cell adhesion, invasion, and migration. TSPAN9 is alluded to have anti-cancer effects in gastric cancer, suppressing proliferation, invasion, and migration in gastric cancer cell lines [57,58]. Despite that, adding TSPAN9 to the COL1A1, COL5A1, ITGA4, and EMILIN1 Cox model further increased the hazard ratio to 36.813 for CAF infiltration in stomach adenocarcinoma samples (Fig. 10a, Table 4). However, the hazard ratios remained zero with the stepwise addition of ITGA4, EMILIN1, and TSPAN9 to the Cox model in adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma ( Table 4). All this data suggested COL1A1, COL5A1, ITGA4, EMILIN1, and TSPAN9 as a poor prognostic CAF signature with high specificity to stomach adenocarcinoma.
We further investigated the expression profiles of the signature genes and their correlation with CAF infiltration (Fig. 10b-e). Mesothelioma and stomach adenocarcinoma displayed a higher expression profile for the COL1A1, COL5A1, ITGA4, and EMILIN1, in comparison to adrenocortical carcinoma and kidney renal papillary cell carcinoma. The expression of TSPAN9 was similar in all four cancers (Fig. 10b). Correlation between the COL1A1 or COL5A1 expressions and CAF infiltration was strong in all four cancers (rho > 0.5). ITGA4 expression poorly correlated with the CAF infiltration, but the correlation was slightly higher in kidney renal papillary cell carcinoma and stomach adenocarcinoma. Strikingly, the correlation of EMILIN1 and TSPAN9 with CAF infiltration was quite strong in stomach adenocarcinoma compared to a poorer correlation in the other three cancers (Fig. 10c).
Further investigation revealed that the expression of ITGA4 increases with stage in stomach adenocarcinoma (Fig. 10d). Although there was a significant decrease in EMILIN1 and TSPAN9 levels in stage 1 compared to healthy stomach tissue, their expression increased again at stage 2, reaching the level of or above that of normal tissues (Fig. 10e-f ). A similar pattern was not observed for adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma (Additional file 2: Fig. S7).
The KM-survival analysis did not indicate a prognostic role for ITGA4, EMILIN1, or TSPAN9 in stomach adenocarcinoma per se (Additional file 2: Fig. S8A-C). However, their hazard ratios increased with the stage (Additional file 2: Fig. S8D-F). This was in parallel to the increase in the poor prognostic impact of CAF infiltration by stage in stomach adenocarcinoma (Additional file 2: Fig. S5), suggesting a stage and CAF dependent role for ITGA4, EMILIN1, and TSPAN9.

Search on drugs that target the poor prognostic CAF signature
Lastly, we searched for currently available drugs that target the five signature genes COL1A1, COL5A1, ITGA4, EMILIN1, and TSPAN9. The drugs that target EMILIN1, and TSPAN9 are not currently available. But our search on DrugBank and DGIB revealed three agents which target COL1A1 and COL5A1: collagenase clostridium histolyticum, halofuginone, and ocriplasmin; and three agents which target ITGA4: natalizumab, firategrast, and BIO-1211 (Fig. 11).
Collagenase clostridium histolyticum and ocriplasmin are enzymes that cleave COL1A1 and COL5A1. They also have proteolytic activity on COL3A1 and FN1, respectively [59,60]. Collagenase clostridium histolyticum is used on skin ulcers to hasten wound healing and Dupuytrens' disease to resolve contractures by digesting collagen [61,62]. Intra-tumoral or intravenous injection of collagenase increased the (See figure on next page.) Fig. 8 Prognostic impact of the six key CAF markers and the CAF infiltration in TCGA cancers. Heatmap of risk scores for A the expression of COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC , B the cancer-associated fibroblast (CAF) infiltration, C the expression of COL1A1 and COL5A1 plus CAF infiltration in different cancers at TCGA dataset. A-C z-score > 0 (p < 0.05) indicates increased risk and z-score < 0 (p < 0.05) indicates decreased risk. ACC: adrenocortical carcinoma, BLCA: bladder urothelial carcinoma, KIRP: kidney renal papillary cell carcinoma, MESO: mesothelioma, and STAD: stomach adenocarcinoma. Data was extracted from TIMER 2.0. D Venn analysis of genes highly correlated (r ≥ 0.5) with COL1A1 and COL5A1 in TCGA samples of ACC, KIRP, and MESO. E The circos plots that show how genes from the input gene lists-25 common genes detected by Venn analysis at (D) (red outer circle) vs. 38 upregulated genes in gastric cancer (blue outer circle)-overlap. The dark orange color at the inner circle represents the genes that appear in both lists and the light orange color represents genes that are unique to that gene list. Purple lines (upper circos plot) link the same genes that are shared by the two lists. Blue lines (lower circos plot) link the different genes which fall into the same ontology term. diffusion of large drug molecules in tumor models [63]. Intraperitoneal administration of collagenase was reported to increase the efficacy of chemotherapy by cleaving the tumor stroma in a rat model of colorectal cancer peritoneal metastasis [64]. Ocriplasmin is used to remove adhesions in symptomatic vitreomacular adhesion [60]. Like our study, another bioinformatics study suggested ocriplasmin as a potential anti-cancer agent [65]. But, to the best of our knowledge, ocriplasmin has not been tested in cancer before.
Halofuginone is an alkaloid that suppresses the expression of the COL1A1 gene, cell migration, and ECM formation. Besides its' antifibrotic and anti-angiogenetic actions, halofuginone shows antiproliferative effects by inhibiting TGFβ/Smad3 signaling [66,67]. Halofuginone, showed an apoptotic effect in prostate cancer and Wilms' tumor cells by inhibiting the transformation of fibroblasts to myofibroblasts [68], which carry similar features to CAFs [7]. Halofuginone also acted synergistically with gemcitabine and suppressed tumorigenesis in a mouse pancreatic cancer model by reducing the number of stromal myofibroblasts and generation of ECM [69].
Integrin α 4 β 1 is a significant therapeutic target in chronic inflammatory diseases and cancer. Natalizumab, the monoclonal antibody against integrin α 4 subtype, was approved in multiple sclerosis and inflammatory bowel disease. However, its long-term use is associated with progressive multifocal leukoencephalopathy [49]. Although abrilumab and vedolizumab are listed as integrin α 4 targeting agents, their action is specific to integrin α 4 β 7 heterodimer [70,71]. Besides monoclonal antibodies, small molecule inhibitors that target integrin α 4 such as firategrast and BIO-1211 are available [72,73]. To the best of our knowledge, these agents have not been tested for their therapeutic efficacy in cancer yet.

Discussion
In this study, we performed a comprehensive bioinformatic analysis to identify poor prognostic CAF markers targeting which may have a therapeutic potential in gastric cancer patients. Our networkbased approach revealed an upregulated ECM protein hub where the CAF markers COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC were the most central genes. High expression of all these CAF markers was associated with high CAF infiltration, tumor progression, mesenchymal phenotype, and decreased survival in gastric cancer patients. Based on the comparative pan-cancer analysis of these key CAF markers and a comprehensive literature search we identified COL1A1, COL5A1, ITGA4, EMILIN1, and TSPAN9 as the signature genes, which potentiated the poor prognostic impact of CAF infiltration specifically in stomach adenocarcinoma. Our findings emphasize the key role of the tumor microenvironment and CAFs in gastric cancer.
Tumor cells dynamically interact with their microenvironment which consists of an ECM compartment and a cellular compartment. CAFs are pivotal cellular components in the tumor microenvironment. They secrete numerous ECM proteins, mainly fibrous collagens (type I, III, and V collagens) and fibronectin. They also remodel the ECM through matrix metalloproteinases (MMPs) which cleave the ECM components, and the lysyl oxidase (LOX) family enzymes which crosslink the collagens. The dynamic remodeling of ECM by CAFs facilitates cancer cell migration and invasion [7,29]. Additionally, CAFs induce ECM stiffness in the tumor microenvironment, which is associated with chemoresistance and poor survival in many cancers [7,29].
CAFs are mostly renowned for their pro-tumorigenic role. However, they can also exhibit an anti-tumorigenic role in a tumor-dependent manner. Whether they exhibit a pro-tumorigenic or an anti-tumorigenic role is highly determined by their secretome [74]. Therefore, identifying the poor prognostic signature of CAFs is of critical importance to develop anti-CAF strategies and Table 3 Parameters of the multivariate Cox proportional regression model for six key CAF markers and CAFs. TCGA data for stomach adenocarcinoma samples were analyzed in the TIMER2.0 immune association -outcome module using the TIDE algorithm for the allocation of samples to the high vs. low CAF infiltration groups. Likelihood ratio and Score log-rank tests were performed (CAF: cancer-associated fibroblast, CI: confidence interval, HR: hazard ratio)  . 9 Enriched ontology clusters for the cancer-associated fibroblast markers and upregulated genes in gastric cancer. A Dendrogram of the enriched ontology clusters (GO/KEGG terms, canonical pathways) for the CAF markers and upregulated genes in GC. B Network representation of enriched terms for a combined list of the CAF markers and upregulated genes in gastric cancer (GC). Each term is represented by a circle node, where its size is proportional to the number of the input genes that fall into that term. The nodes are presented as pie sectors where each pie sector is proportional to the number of hits originating from the CAF markers list (red) and upregulated genes in GC (blue). C The gene ontology terms for the same network nodes in B where each color represents different cluster identities given in the label. D Representation of the same network nodes in B and C colored by p-value, as shown in the color scale. B-D The black arrows show the nodes that fall into the Integrin α4β1 pathway in GC upregulated genes list. Only one node that falls into the Integrin α4β1 pathway shown with the red arrow was common in both the CAF markers list and the GC upregulated genes list. The data was analyzed, and the network layouts were prepared in Metas cape. org select the patient populations that will benefit from these modalities. Our study and others' findings demonstrated a poor prognostic role for CAFs in gastric cancer [8,9] (Table 3 and Additional file 2: Fig. S5). Zeng et al. computationally analyzed the cell infiltration pattern in the tumor microenvironment of 1524 gastric cancer patients. The fibroblast infiltration was the greatest risk factor in tumor microenvironment phenotype with the poorest overall survival [10]. Despite that, the poor prognostic signature for CAFs is not clear in gastric cancer, and drugs that target CAF-specific protumorigenic processes are not available in the clinic. Hence, in this study, we aimed to address these points.
To identify a poor prognostic signature for CAF infiltration in gastric cancer, we first identified the key CAF markers with a network-based approach in gastric cancer. We identified ECM components COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC as the pivotal CAF markers in gastric cancer, which were separately reported as biomarkers of gastric cancer in different studies [75][76][77][78]. We further validated their differential expression and poor prognostic significance in gastric cancer by analyzing the TCGA, GTEx, and ACRG cohorts. Additionally, we demonstrated that the high expression of these six CAF markers is associated with mesenchymal phenotype gastric cancer, which has a poorer prognosis and response to chemotherapy than the gastric tumors with the epithelial phenotype (Fig. 6). The higher content of stroma in mesenchymal phenotype gastric cancers and the role of EMT in CAF generation may explain this association. Interestingly, we could not observe a similar association in the diffuse histopathological subtype of gastric cancer which is more related to a mesenchymal phenotype compared with the intestinal histopathological type. However, Oh. Et al. suggested that the mesenchymal phenotype gastric cancers are a subset of diffuse gastric cancers [25]. Hence the heterogeneity of diffuse gastric cancers in terms of a mesenchymal or epithelial phenotype may explain this discrepancy.
In the next step, we performed pan-cancer profiling of the six key CAF markers. Although we identified five different cancers in which all the six CAF markers were upregulated, these markers or CAF infiltration did not exhibit poor prognostic effect in these cancers. These observations suggested that the six CAF markers we identified may have more critical roles in the poor prognostic impact of CAFs in gastric cancer. Then, Table 4 Parameters of the stepwise multivariate Cox proportional regression model for cancer-associated fibroblasts and signature genes. TCGA data for stomach adenocarcinoma, adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma were analyzed in the TIMER2.0 immune association -outcome module using the TIDE algorithm for the allocation of samples to the high vs. low CAF infiltration groups. Likelihood ratio and Score log-rank tests were performed (CAF: cancer-associated fibroblast, CI: confidence interval, HR: hazard ratio) we searched for other cancers in which all the six CAF markers and CAF infiltration are poor prognostic. We established a poor prognostic gene signature based on the comparative analysis of gastric cancer with these tumors, namely adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma.
Among the six CAF markers, we identified COL1A1 and COL5A1 as the two markers which potentiated the poor prognostic impact of CAF infiltration most in gastric cancer ( Table 3). The opposite effect of this dual gene signature in adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma (Fig. 8e) led us to identify ITGA4, EMILIN11, and TSPAN9 as interactors that Currently available drugs that target COL1A1, COL5A1, and ITGA4. Drugs (shown with blue squares) listed in DrugBank and Drug-gene interaction database that target COL1A1, COL5A1, or ITGA4 (shown with orange circles). Among the ITGA4 targeting drugs, the action of abrilumab, and vedolizumab, which are shown in light blue, are specific to integrin α4β7 rather than integrin α4β1. Natalizumab, firategrast, and BIO-1211 target integrin α4 potentiate the poor prognostic effect of CAFs specifically in gastric cancer. To the best of our knowledge, this is the first time the ITGA4, EMILIN1, and TSPAN9 are put forth as poor prognostic signature genes for CAF infiltration in gastric cancer.
Recently, Liu et al. suggested TGFB2, VEGFB, COL10A1, AREG, and EFNA5; and Grunberg et al. suggested THBS1, THBS2, and INHBA as poor prognostic signatures for CAF infiltration in gastric cancer [8,32]. But the z-score and hazard ratio for the dual COL1A1 and COL5A1 gene signature we identified was higher compared to both signatures (Table 4, Additional file 1: Fig. S8). The hazard ratio of CAF infiltration in COL1A1, COL5A1, ITGA4, EMILIN1, and TSPAN9 signature was even higher (Table 4), presenting a 36.8 times higher risk of death in gastric tumors with high CAF infiltration. One possible advantage of THBS1, THBS2, and INHBA signature maybe its predictive ability in liquid biopsies since the signature genes are secreted through extracellular vesicles [32]. Further studies are needed to assess and compare the predictive potential of all these signatures in different biopsy specimens.
The delineation of ITGA4 as a poor prognostic factor was not surprising since increased expression of integrin α 4 β 1 is associated with tumor progression and ECM components secreted by CAFs bind integrins to activate pro-tumorigenic pathways [41,79,80]. However, the poor prognostic effect of EMILIN11 and TSPAN9 was surprising since EMILIN-1 is regarded as a tumor suppressor which acts synergistically with TSPAN9. Knockout of EMILIN11 or suppression of EMILIN-1 -integrin α 4 β 1 interaction was associated with decreased expression of the tumor suppressor PTEN, and increased activity of PI3K/Akt and ERK1/2 pathways, leading to hyperproliferation of dermal fibroblasts and keratinocytes, increased skin carcinogenesis and lymph node metastasis [50,53]. Knocking out EMILIN1 or transgenic expression of an EMILIN1 mutant with impaired binding to integrin α 4 β 1 increased the susceptibility of mice to develop colon cancer [81]. Based on these findings EMILIN1 was proposed as a tumor suppressor. However, EMILIN1 overexpression was detected in serous ovarian carcinoma, soft tissue osteosarcoma, and low-grade glioma (LGG) which are malignant tumors with high recurrence rates [54,55,82]. This suggests a tissue-dependent anti-tumorigenic or pro-tumorigenic role for EMILIN-1.
Recently, EMILIN-1 was suggested to increase TSPAN9 expression in gastric cancer cell lines and form a complex with TSPAN9 to synergistically inhibit FAK/Ras/Erk pathway and suppress invasion and migration. Since overexpression of only EMILIN1 did not induce a similar anti-tumor response in gastric cancer cell lines, it was suggested that the anti-tumor effect of EMILIN-1 may be dependent on TSPAN9 [56]. Therefore, we tested whether the addition of TSPAN9 to the CO1A1, COL5A1, ITGA4, and EMILIN1 Cox model can decrease the poor prognostic impact of CAF infiltration in gastric cancer. Unexpectedly, the hazard ratio and risk score for the CAF infiltration increased further in stomach adenocarcinoma but remained zero in adrenocortical carcinoma, kidney renal papillary cell carcinoma, and mesothelioma (Table 4).
Although some studies suggest an anti-cancer role for TSPAN9 [57,58], high TSPAN9 expression was associated with resistance to 5-fluorouracil via suppressing autophagy in gastric cancer cell lines [83]. Expression of TSPAN9 was significantly lower in gastric cancer tissue compared to adjacent normal gastric tissue in a cohort of 105 gastric cancer samples. However, the same cohort reported high TSPAN9 expression as a poor prognostic factor for survival [84]. Therefore, whether EMILIN1 and TSPAN9 exert an anti-tumor effect or protumorigenic effect in gastric cancer is not clear yet.
Our investigation of TCGA stomach adenocarcinoma data does not suggest either a poor or a good prognostic role for EMILIN1 or TSPAN9 per se (Additional file 2: Fig. S8B-C). Despite that, the EMILIN1 and TSPAN9 displayed a stage-dependent increase in expression (Fig. 8f ), and their hazard ratios increased by stage in parallel to the stage-dependent increase in the hazard ratio of CAFs (Additional file 2: Fig. S8E-F). Moreover, the stepwise addition of these two genes as covariates to the CO1A1, COL5A1, and ITGA4 multivariate Cox model tremendously increased the hazard ratio for the poor prognostic impact of CAF infiltration, which suggested a stage and CAF dependent role for EMILIN1 and TSPAN9 in gastric cancer (Table 4). This may explain why these two genes display anti-tumor effects in monoculture gastric cancer cell lines but display poor prognostic effects in the gastric cancer patient cohort by Feng et al. [84] and TCGA stomach adenocarcinoma cohort we analyzed in this study.
It may be speculated that the ECM remodeling enzymes secreted from CAFs may cleave EMILIN-1 and prevent its anti-tumor action together with TSPAN9 despite their high expression in the tumor tissue. Accordingly, proteolytic cleavage of EMILIN-1 by MMPs or neutrophil elastase was suggested as a mechanism for pro-tumorigenic effect in some tumors with high EMILIN1 expression [48,85,86]. For ECM proteins, it is also not unusual to serve different functions in cleaved forms vs. multimeric forms [87]. A similar mechanism may explain the context-dependent role of EMILIN-1 and TSPAN9 in gastric cancer. Although there is no evidence in the literature yet for multimerization of EMILIN-1, this may also be a possible mechanism for its context-dependent action, since our analysis pointed out "protein-crosslinking" as an enriched biological process in gastric cancer. Elucidation of these molecular and cellular mechanisms may present EMILIN-1 and TSPAN9 as new therapeutic targets in gastric cancer. This will be addressed in our future studies.
To build the poor prognostic gene signature for CAF infiltration, we analyzed the TCGA stomach adenocarcinoma data in TIMER 2.0. TIMER 2.0 implements four different algorithms, namely EPIC, MCP-Counter (Microenvironment Cell Populationscounter), Xcell, and TIDE (Tumor Immune Dysfunction and Exclusion) to predict the relative proportions of different cell populations in tumor samples [88][89][90][91]. All these algorithms device reference gene expression profiles for each cell type, established from the RNA-seq profiles of circulating immune cells and non-cancerous cells that infiltrate the tumors. One limitation to the study may be that only the TIDE algorithm predicted a poor prognostic impact for CAF infiltration in COL1A1, COL5A1, ITGA4, EMILIN1, and TSPAN9 multivariate Cox model in gastric cancer. One reason may be that the reference gene expression profiles for CAFs in all four algorithms are different, which leads to differences in the allocation of samples to high vs. low CAF infiltration groups. Dissecting these differences in detail may improve the power of these algorithms to predict the abundance of CAFs in tumor samples. Moreover, mostly melanoma samples are used to establish reference gene expression profiles for tumor-infiltrating cells. These reference gene expression profiles may not exactly reflect the gene expression pattern for CAFs in gastric cancer or other cancers. Besides intertumoral heterogeneity, intra-tumoral heterogeneity in CAFs further complicates the picture [7]. Building tumor and subclone-specific reference gene expression profiles for CAFs may better illuminate the prognostic role of CAFs and their interactor genes in cancer. Such an approach may also reveal new molecular targets to prevent CAF infiltration in cancer.
Among the signature genes we identified in this study, targeting the COL1A1, COL5A1 and ITGA4 may have a therapeutic potential in gastric tumors with high CAF infiltration. Our search in drug databases brings out collagenase clostridium histolyticum, halofuginone, and ocriplasmin as agents that act on COL1A1 and COL5A1. Their anti-cancer action should be validated first in gastric cancer cell models and in vivo studies. Among these three, halofuginone seems to be the closest candidate for use as an anti-cancer agent in gastric cancer since it showed promising anti-cancer effects in other cancers. It may also prevent the formation of CAFs [68,69]. Moreover, collagenase clostridium histolyticum and ocriplasmin carry risks, since cleavage of collagens may lead to a release of several growth factors that induce tumor progression. Additionally, their systemic use can be problematic in cancer due to the risk of organ and vascular toxicity [63]. Active targeting of the gastric tumors via nanocarriers or local administration may allow their use in gastric cancer. But still, there may be destructive effects on healthy gastric tissue limiting the doses that could be used in patients safely. All these points should be addressed in future studies.
Integrin α 4 β 1 also has potential as a therapeutic target in cancer since it has a key role in metastasis, chemoresistance, angiogenesis, and lymphangiogenesis [41]. Despite that, agents targeting integrin α 4 β 1 have not entered the clinical trials for cancer. The concerns about the risk of cancer development with targeting integrin α 4 β 1 may be the reason since this action inhibits the migration of lymphocytes. However, the comparative analysis did not reveal an increased risk of cancer with natalizumab [92]. Uncovering the signaling mechanisms by which integrin α 4 β 1 contributes to cancer may lead to the development of new strategies for targeting integrin α 4 β 1 without raising the concerns about cancer development.

Conclusion
In this study, we identified the six key CAF markers namely COL1A1, COL1A2, COL3A1, COL5A1, FN1, and SPARC in gastric cancer that can be used as predictors of CAF infiltration and prognostic biomarkers in gastric cancer. With further analysis, we revealed COL1A1, COL5A1, ITGA4, EMILIN1, and TSPAN9 as a poor prognostic gene signature for CAF infiltration, with high specificity to stomach adenocarcinoma. This signature could be translated to the clinic with further studies as a predictive tool for poor prognosis. Testing the candidate drugs, we identified in this study, with further in vitro and in vivo studies may present them as potential drugs for the treatment of gastric cancer. More importantly, investigating the mechanisms by which the signature genes strengthen the poor prognostic impact of CAFs may put forth new molecular targets for the effective treatment of gastric cancer. These points will be addressed in our future studies.